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Abstract. We investigate the dynamics of a single-ended TV-state molecular zipper 
^ . based on a model originally proposed by Kittel. The molecule is driven unidirectionally 

t^ I towards the completely unzipped state with increasing time t. The driving lowers the 

energies of states with k unzipped links by an amount proportional to kt. We solve the 
Pauli rate equation for the state probabilities and the partial differential equations, 
which yield the probability distributions for the work performed on the zipper and for 
' y i the heat exchanged with the thermal reservoir. Similarly to the related equilibrium 

model, two different regimes can be identified at a given temperature with respect to 
(^ ' released molecular degrees of freedom per broken bond. In these two regimes the time 

evolution of the state probabilities as well as of the work and heat distributions show 
^S| ' a qualitatively different behavior. 
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1. Introduction 

The Watson-Crick double-stranded form represents ttie tlier mo dynamically stable state 
of DNA in a wide range of temperature and salt conditions. However, even at standard 
physiologic conditions, there always exists possibility that the double-helix is locally 
unzipped into two strands both at its ends and in its interior [H El El HI [5] . If the interior 
unzipping is neglected, the unfolding of the two strands can be described by a simple 
zipper model [HI [71 [8] . This model, while including several simplifications, emphasizes 
the essential ingredient of the unzipping process, that means the competition between 
the entropic forces which tend to open the macromolecule, and the energetic forces that 
tend to condense it into its double-stranded form. 

In the last two decades, new experimental techniques have been developed for 
detailed analysis of unzipping processes. A powerful technique are single molecule 
manipulations by optical tweezers [9l [TOl HH [121 HSl HI] ) where biomolecules are unfolded 
and refolded by applying mechanical forces. Because single molecules are subject to 
strong fluctuations, a stochastic description of the unfolding/refolding kinetics becomes 
necessary. Investigation of these stochastic processes can be useful to understand how 
biomolecules unfold and fold under locally applied forces [ISl [161 [TTl [E] , ^-s, for example, 
when mRNA passes through the ribosome during the translation process [191 120], or 
when DNA is unzipped by helicase during the replication process [211 [22] . 

A particularly interesting part in the analysis of unzipping processes is the 
application of integral or detailed fluctuation theorems [231 1211 [251 El [26] to estimate 
free energy differences between folded and unfolded states (see, e.g., [2Z])- Their 
advantage is that they can be applied also to protocols, which drive the considered 
system far from equilibrium. It is thus not necessary to perform the unfolding/folding 
under quasi-static near-equilibrium conditions, where the process becomes reversible. 
Most popular theorems are the Crooks fluctuation theorem [28] and Jarzynski equality 
|29j . These relate to the distribution of work performed on the molecule during the 
process. Unfortunately it is not easy to get theoretical insight into characteristics of 
the underlying work distributions in far-from-equilibrium processes, since the work is 
a functional of the whole stochastic trajectory. Investigations have been conducted for 
simple spin systems driven by a time- dependent external field [301 SB ESI EHl [SI [35] and 
for diffusion processes in the presence of a time- dependent potential [361 [37]. Analytical 
solutions are known for two-state systems [32l [331 EH [35] and for systems with one 
continuous state variable |36l [37]. 

In this study we present analytical results for the work distribution for a multi-state 
system, which is motivated by a model proposed by Kittel for describing the melting 
transition of DNA molecules [8j. In order to derive analytical results, we had to consider 
a stochastic process with directed forwarded and forbidden backward transitions between 
states. As a consequence, detailed balance is broken and fluctuation theorems, as 
mentioned above, will not hold true. On the other hand, kinetic Monte Carlo simulations 
of the stochastic process show that, for at least certain parameter settings in the model. 



the restriction of forbidden backward transitions is not so severe and events dominating 
the integrand in the Jarzynski equahty are not very rare. It is important to point out 
that these parameter settings are not related to any experimental conditions. In more 
realistic settings rare events would play the decisive role and in such cases it becomes 
difficult to determine tails of the work distribution with sufficient accuracy. Also the 
relation of the work considered in our study to the thermodynamic work measured in an 
unzipping experiments needs to be treated with care. These problems are discussed in 
detail in Sec. [2] and they imply that the theory cannot be applied to experiments at this 
stage. Our findings should nevertheless be useful because connection to experiments 
seems not completely out of reach and because they widen the range of hopping models, 
where analytical results for work distributions are available. 

2. Unzipping in an extended Kittel model 

Kittel's model [8] is a simplified version of the Poland-Scheraga model [1] for the 
equilibrium properties of DNA molecules, which got renewed attention in the last ten 
years [39]. In contrast to the Poland-Scheraga model, it disregards the possibility of 
any interior openings (bubbles) of double-stranded parts of the molecule. Despite this 
simplification, it is already sufficient to understand the origin of a melting transition. 

A molecule in Kittel's model [8] has N links in its fully folded state, see Fig. [TJ 
Different states k = 1,...,N oi the molecule refer to configurations, where {k — 1) 
of the links are opened in a row. The difference Ek — Ei of the internal energy of a 
state with {k — 1) open links and the ground state k = 1 (no open hnks) is equal to 
{k — 1)A, corresponding to a loss of chemical bond energy A per broken link. With 
each broken link, the molecule gains G degrees of freedom, which characterize the win of 
conformational degrees of freedom when double-stranded DNA is transferred to single 
stranded DNA. The entropy difference Sk — Si between the state k and the ground state 
then becomes /cbIhG'^"^ = kB{k — l) InG, where k-^ is the Boltzmann constant. The free 
energy F^ of state fc at a temperature T is thus given by Fk = Fq + (A; — 1) [A — k^T In G] 
and the equilibrium properties can be readily worked out by considering the partition 

sum Z = Z]fc=i6^P(~/^-^fc)' where /3 = l/k^T. 

In extending this model to treat unfolding kinetics, the molecule is supposed to 
unzip successively, one link in each step, as a consequence of a strong external driving. 
Such external driving can be caused by a local force or a pH gradient, as indicated in 
FlG.lll With respect to pulling experiments on single DNA molecules out of equilibrium, 
there is evidence that a neglect of interior openings can become even less relevant than 



Figure 1. The zipper model. 



for the equilibrium properties. For example, far away from the melting temperature, 
e.g. at standard room temperature conditions (298 K), the unfolding kinetics of DNA 
hairpin molecules could be successfully described by assuming no interior openings and 
thermally activated unzipping transitions [161 [38] • Moreover, molecular fraying can be 
identified in individual unzipping trajectories by considering the size of force jumps. 
Thereby states with interior openings can be systematically excluded from the analysis 

m- 

We assume that the driving lowers the energy differences by an amount proportional 
to {k — l)t, i.e. open states with a larger number of broken links are favored with 
increasing time. If we take the ground state energy as the reference point, Ei = 0, we 
obtain for the free energies -Ffc(t) = -Efc(^) ~ kBTlnG^^^ of the states 

Fk{t) = {k-l)[A-vt-kBT\nG], k = l,...,N, (1) 

where the parameter v has the dimension of an energy rate and characterizes different 
speeds of the unfolding in response to different strengths of the external driving. 

Following established theoretical descriptions for the rupture kinetics of the bonds 
[38| HHl E] , the Kramers-Bell form |12] is used for the time- dependent transition rate 
^k,k+i from state k to state [k + 1), 

Xk,k+iit) = z>exp[-/3Ffe,,,+i(t)] . (2) 

Here 9 is an attempt frequency and -FA;,fe+i(^) denotes the free energy barrier for the 
transition, i.e. the difference of the free energy at the saddle point separating states 
k, {k + 1) and the free energy Fk in state k. The barrier Fk^k+iif) is considered to 
be composed of a bare, fc-independent free energy barrier, F^, which is modified by an 
amount proportional to the free energy difference, [Fk+i{t) — -Ffc(t)], 

Fk,k+i{t) = F, + ^[Fk+i{t) - Fk{t)] = Ffe + 7[A - t;t - k^TlnG] . (3) 

In the following we set 7 = 1 [43j. Note that F^^k^iit) and thus \(t) = \k,k+i{t) in 
Eq. dJ]) are independent of k. The attempt frequency z> in Eq. (jj]) was reported |11] to 
be approximately proportional to the ratio of the diffusion constant of the molecule in 
water to the water viscosity. We assume here a linear dependence of this ratio on the 
temperature T, i.e. we take 9 = z/(T/To), where u and Tq are positive constants. Thus 
we obtain 

A(t) = Xk,k+iit) = gexp[-P{d - vt)] , (4) 

where 

g = uG^, d = A + F,. (5) 

-'o 

If detailed balance is obeyed, i.e. Xk,k+i(t) Gxp[—/3Fk(t)] = Afc+i,fc(t) exp[— /3Ffc+i(t)], 

the rates \k+i,k(t) = Xt for backward (refolding) transitions become both independent 

of k and t, 

Xb = Xk+i,k{t) = —exp{-f3Fb) . (6) 



Note that G appears here in the denominator, which means that the ratio Xb/X{t) = 
exp(/9A — (5vt)/G of backward to forward rates becomes small for large degeneracy 
factors G. This reflects the fact that it is difficult for the flexible unfolded part of 
the molecule to find proper configurations, which would allow for a reformation of 
(hydrogen) bonds. 

If the model would refer to the hopping motion of a particle between time-dependent 
energy levels Ek{t), the work performed on the system (for one realization of the 
stochastic process) would be 

w(t)= / dt' VEfc(05fcD(t') 

•^0 k=l 

D(t)-1 

= [^D(t)(t) - ^D{t)(tD{t))] + Yl [^'^('^'i+l) - M'tk)] 

k=l 
D(t)-1 

= [FD(t)(t) - FD(t)(tD(t))] + J2 l^tk+i) - Fkih)] , (7) 

fe=l 

where D(t) denotes the (random) state of the molecule at time t, Ek(t) = dEk(t)/dt, and 
tk is the (random) time at which the transition from state k to {k + 1) (rupture of fcth 
bond) takes place (ti = being the initial time). The last line in Eq. ([7]) follows from 
the fact that the entropy in the extended Kittel model is not dependent on time, i.e., 
for given k, differences between internal energies and between free energies at distinct 
times are equal. 

The work in Eq. ([7]) can be related to the thermodynamic work Wj in an unzipping 
experiment under force control. As pointed out in ^5], work in thermodynamics is 
the internal energy transferred to a system upon changing the control parameters for 
given system configuration. This means that, when the force / = f{t) is the control 
variable and the molecular extension x the conjugate configurational variable, one has 
Wf{t) = — Jr,J dfx |16]. Stochastic changes of the molecular extension occur mainly 
due to bond rupture, while in between transitions the response will be rather smooth 
and, under neglect of small thermal fluctuations, can be represented by a function 
x{k, /). This specifies the mean end-to-end distance of the unfolded part in state k at 
force / (as often modeled by, e.g., the freely jointed or worm-like chain models from 
polymer physics). With /^ = /(t^) we then have 

Wf{t) = - dfx = ~ dfxiDit)J)- J2 dfxikj). (8) 

Jf{0) Jfoit) k=l -^Sk 

In more detailed energy landscape models, the free energy Fioi{k, /) of the molecule 
in state k under loading / can be represented as -Ftot(fc, /) = FQ^k) — fx{k, f) + Fstr{k, f) 
(see, e.g., [E]), where Foi^k) is the free energy in the absence of loading and Fstr{k, f) 
the elastic energy of the unfolded part upon stretching. For smooth response under 



stretching the latter is given by 

Fstrik,f)= r ' dxfik,x) = fxikj)- [ df'xikj'), (9) 

Jx{k,0) Jo 

where f{k,x) denotes the inverse function of x{k, f) with respect to /. Accordingly, at 
a given k, the work for stretching upon increasing the force from fa to ff, becomes 
rfb 

- / dfx{k, f) = F^tr{k, fb) - fbX{k, ft,) - F^tr{k, fa) + faX{k, fa) 

Jfa 

= FtotikJb) - FtotikJa) . (10) 

This just means that the stretching at fixed k is assumed to take place quasi-statically, 
i.e. the variation of the force is supposed to occur on a time scale much slower than the 
correlation time of end-to-end distance fluctuations. Inserting this result into Eq. (IHl) 
yields 

D(t)-1 

Wf{t) = -[FUk, fit)) - F,^,{k, D(t))] - J2 l^^otik, /(tfc+i)) - Ftot(fc, fih))] , (11) 

fc=i 

from which it becomes clear that Wf equals {—W), if Ftot{k, f{t)) is identified with the 

free energy Fk{t) in the extended Kittel model. 

As said in the Introduction, we were able to find analytical results for the work 
distributions, if the backward rates in Eq. ([6]) were negligible. We are interested here 
in the model itself and will not make attempts to assign values to the parameters in 
the transition rates A(t) [Eq. (jl])] and Af, [Eq. (Q], and appearing also in the state 
free energies in Eq. ([T]), which are connected to real experiments. Nevertheless, with 
respect to the value of the analytical results in connection with fluctuation theorems, 
the question arises whether a neglect of backward rates could be acceptable, at least 
for certain parameter settings. To check this, we have performed kinetic Monte Carlo 
simulations of the stochastic process. Because the forward rate from Eq. (jl]) can be 
written as \(t) = XbGexp{—f3A + f3vt) and the free energy differences appearing in 
Eq. ([7]) by [-Ffc(tfc+i) — Ffc(tfc)] = —{k — l)v(tk+i — t^), the dynamics and energetics 
of the model is completely specified by the parameters Xb, G, v, and T (and A^ if we 
consider a complete unfolding). For illustration and discussion of representative results, 
we here use A, A/ks, and u'^ as units for energy, temperature and time, respectively. 

Simulations were performed for fixed A^ = 10, T = 1 and A;, = 1, and a set of v and 
G values varying in the intervals v = 0.01 — 3.3 and G = 10 — 1000, respectively. 
We always started the unfolding from the fully closed state, i.e. Pk{0) = 6k,i |17] . 
Probabilities PN(t) of complete unfolding (occupation of state A^) until time t were 
determined and an unfolding time tu defined by requiring that at t = tjj the zipper 
has unfolded with a probability of 99.9%, i.e. Pwitu) = 0.999 (see Sec. |5]). We then 
considered the work distributions p{w,tu) at t = tjj and the weighted distributions 
exp(— /3w)p(w, t(/), corresponding to the integrand in the average (exp(— /Sw)) = 
Jdw exp(— /3w7)p(w, tf/) as it appears in the Jarzynski equality. It was found that 
backward rates turn out to have a minor importance when the forward rates are much 
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Figure 2. Simulated results of the model presented in Sec. [2] with backward transition 
rate Xt — 0.133 (blue lines) in comparison with analytical results (see Sec. HJ when 
neglecting backward transitions (green circles). Panels labeled with a) and b) refer 
to degeneracy factors G = 10 and G — 1000, respectively. The remaining parameters 
are T = 1, v = 0.25, and N — 10. Panels al) and bl) show the probability PN{t) 
that the zipper has unfolded up to time t and the insets depict the time-dependence 
of the forward transition rates (full line); the constant backward rate is indicated by 
the dashed line. Panels a2) and b2) display the probability densities p{w,tij) at the 
unfolding time tjj {tjj — 10.67 for G = 10 and tjj — 0.47 for G = 1000), and panels a3) 
and b3) the weighted probability densities exp(— /3w)p(w,i[/). 



larger than the backward rates during the whole unfolding process or when v is large 
enough. Specifically, the error in calculating (exp(— /3w)) is smaller than 5% when 
t; > 3 for G = 10, t; > 0.1 for G = 100, and v > 0.01 for G = 1000. As representative 
examples we show in Fig. [2] simulated results (blue lines) with nonzero backward rates 
in comparison with analytical results (green circles) for zero backward rates (see Sec. H]) 
for PAr(t), p{w,tu), and exp{—(3w)p{w,tu), and parameters v = 0.25, G = 10 [panels 
labeled with a)] and G = 1000 [panels labeled with b)]. As can be seen from the 
figure, for the case of large G (small backward transitions), the simulated results for 
exp{—f3w)p{w,tu) are almost indistinguishable from the analytical results. 

The remaining part of the paper is organized as follows. In Sec. |3] we specify 
the equations for the time evolution of the state probabilities [Eq. f fT2|) ] and for the 
quantities describing the energy transformations [Eq. (TT71)]. In Sec. IHwe derive exact 
analytical solutions of these equations and in Sec. |5]we discuss our findings. 



3. Time evolution of state probabilities and work 



L(t) 



v 



(13) 



/ 



Let Pkit), k = 1,. . . ,N, be the occupation probabilities of the k-th state. The time 
evolution of these functions is governed by the Pauli rate equation with transition rates 
given by Eq. (|4]). Formally speaking, the unzipping process is described by the time- 
inhomogeneous Markov process D(t), where D(t) = k if the system resides in state k at 
time t. The Pauli rate equation can be written as 

— M(t) = hit) R(t) , M(0) = I , (12) 

at 

where I is the (A^ x A^) unity matrix, L(t) is the (A^ x A^) matrix of the transition rates, 

/ -A(t) \ 

A(t) -A(t) ••• ; 

••• ■•• ■•• ; 
; ■•• ■•• -A(t) 

... A(t) 

is the (A^ X A^) matrix of the transition probabilities with the matrix elements 

Rij{t) = Prob { D(t) = 1 1 D(0) = j } . (14) 

The matrix of transition probabilities evolves an arbitrary column vector of the initial 
occupation probabilities, p(t) = M(t)p(0). In the following, we always start with 
the completely closed zipper, that means Pfc(O) = Ski- The individual occupation 
probabilities then are 

Pi{t)=Ra{t). (15) 

According to Eq. (jTj) the work W(t) is a functional of the process D(t). For 
an analytical treatment it is useful to introduce the augmented process {W(t), D(t)} 
[i9l [50l [33] which describes both the work and the state variable. This augmented 
process is again a time non-homogeneous Markov process and its one-time properties 
are described by (A^ x A^) matrix G{w,t) with the matrix elements 

Prob { W(t) e{w,w + e) and D(t) = i \ W(0) = and D(0) = j } 



and 



Gij{w,t) 



lim- 



The time evolution of G{w, t) is given by [49l |50l [33] 



|c(.,t) 



-AE(t)+L(t) 



G{w,t) 



(16) 



G{w,0) = 5{w)I. (17) 



m 



Here E(t) is the diagonal matrix E(t) = diag{Ei(t), . . . , EN{t)}. Notice that 
J^^dwG{w,t). Eq. f[T7j) represents a hyperbolic system of A^^ coupled partial 
differential equations with time-dependent coefficients. Its exact solution will be given 
in the following Sec. [H 

The matrix G{w,t) provides a complete description of the energetics of the 
unzipping process. The joint probability density for the internal energy U(t) = Eo{t){t) 



9 

and the work W(t) performed on the system during the time interval [0, t] (regardless 
of the final state of the system at the time t) is given by 

N 

^u, w,t) = J2^ b^ - Mt)] Gii{w, t) , (18) 

where S{x) is the Dirac 5- function. The last function already yields the probability 

density for the work, 

/OD 
du ^iu,w,t) . (19) 

-oo 

An analogous integration over the work variable w gives the probability density of U(t). 
Furthermore, the first law of thermodynamics implies U(t) — U(0) = U(t) = W(t) + Q(t) 
(note that U(0) = for our setting) and, accordingly, ^{u, w, t) gives also the probability 
density of the heat Q(t) transferred from the reservoir during the time interval [0, t] |35j . 
The mean values of the internal energy, work and heat, are then 

U{t) = / dudwu^{u,w,t) = YEiit)piit) , (20) 



oo ^oo 



/oo /"OO ^^ 

/ dudwu^{u,w,t) = y ^Ei{t)pi{t) 
-OO J —oo -1 

dudww ^{u,w,t) = y^ / dt Ei{t)pi{t) 

/oo ^OO -" ft 

/ dudw {u — w) ^{u,w,t) = y^ / dt Ei{ 
OO J —oo -1 Jo 



W{t)= / / dudww^{u,w,t) = y / dtEi{t)pi{t) , (21) 



— oo J — oo 



Q{t) = dudw{u-w)^{u,w,t) = }_^ dt Ei{t)pi{t) . (22) 

J —oo J —oo -1 Jo 

Notice that these mean values can be also calculated directly from the solution of the 
Pauli rate equation (TT2l) . However, for higher moments, we already need the function 
C,{u,w,t). For example, the variances discussed in Sec. [5] are given by 

/oo /»oo 
/ du d^ u^ ^{u, w, t) - [f/(t)]2 , (23) 

oo J —oo 

/OO /'OO 

/ du dw w^ ^(m, w, t) - [W{t)f , (24) 

OO J —oo 

/•oo /*oo 

[Ag(t)]2 = / / dudw{u-w)^^{u,w,t)-[Q{t)]\ (25) 



oo ■/ — oo 



4. Solution of the model 

Owing to the simple two-diagonal structure of the matrix L(t) in Eq. ( 1T3l) . the system 
( 1T2|) can be solved by simple integrations. By defining 



/ dt"X{t") = -^ [exp{/3vt) - exp{(3vt')] 
Jt' P'^ 



A{t,t')= dt"X{t") = — [exp{(3vt)-exp{Pvt')], (26) 
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where a = gexp{—l3d), a recursive treatment of Eq. ( 1121) yields 

%(t) = exp[-A(t,0)] , j = l,...,N , 

.t (27) 

Ri-{t) = / dt' exp [-A(t, t')] A(t') Ri-ij{t') , j <i< N , 

Jo 

RNjit) = [ dt'A(t') RN^ijit') , j = l,...,N . 

Jo 
When solving these recursive relations, we obtain the lower triangular matrix with the 

nonzero matrix elements 

Rj+kjit)=^-^^^^^exp[-Ait,0)] , j = l,...,iV-l; A; = 0,...,iV-l-j, (28) 

[Aft 0)1 ^"-^ 
^^^•W = (V ^jy ^^P i-Ht, 0)] ii"i(l, iV + 1 - j; A(t, 0)) , j = l,...,Ar. (29) 

where iFi{a,b]x) denotes the confluent hypergeometric function. 

For solving Eq. ( TT71) we first perform a Laplace transform with respect to the 
work variable w. To keep the notation simple, we use the same symbols for the original 
functions and for the transformed ones. The transformed functions will be distinguished 
by explicitly giving the complex variable s conjugate to w. After performing the Laplace 
transformation, we get the system of ordinary differential equations 

d r • 1 

—G{s,t)= -sE(t) + L(t) G{s,t) , G(s,0)=I. (30) 

at I J 

The matrix which multiplies G(s, t) on the right hand side is again a lower two-diagonal 
one. Therefore, similarly to Eq. fl27|) . we find the recursive relation 



Gij{s,t) = , i <j , 

G,,{s, t) = exp [-A(t, 0)] exp {-s[E,{t) - E,{0)]} , j = l,...,N 

Gij{s,t) = [ dt'exp [-A{t,t')] exp {-s[Ei{t) - E,{t')]} A(t') G,-y{s,t') , j<i<N , 

Jo 

GNjis,t)= [ dt'exp{-s[Er,{t)-EN{t')]}X{t')Gr,-i,is,t') , j = l,...,Ar 

Jo 

Notice that the matrix G{s, t) is again a lower triangular one. We now want to solve 
these recursive relations. It turns out that all matrix elements of the matrix G{s,t), 
except the matrix elements Gj\fj{w, s), j = 1, . . . ,N — 1, can be explicitly evaluated by 
simple integrations: 

\ ,n_ ^ [exp(/3t;t) - exp(si;t)] \ 
G,+.,(s, t) = l^P-^ _ L- exp [s{j - l)vt] exp [-A(t, 0)] , (32) 



(31) 



11 

for j = N, k = and also for j = 1, . . . ,N — 1, k = 0, . . . ,N — 1 — j. Moreover, we were 
able to carry out the inverse Laplace transformation of these functions. The resulting 
diagonal elements Gjj{w,t), j = 1,. . . ,N are proportional to Dirac (5-functions. The 
remaining ones possess a finite support, i.e. they are proportional to the differences of 
the unit-step functions 0(a, b;x) = Q{x — a) — Q{x — b) , 



Gj,{w, t) = 6[w + (j - l)vt] exp [-A(t, 0)] , j = 1, . . . , iV 



Gj+kj{w,t) = exp[-A(t 



a 



V 



exp {f3[w + {k + j - l)vt]} X 



(fc-1)! 
^ y (-1)' Q[il-J-l)vt,il-j)vt;w] 



ii-fc 



(33) 



(34) 



^/!(fc-/)! [^ + (^- + /_i)^t]^ 

These expressions are valid for j = 1 , . . . , A^ — 1 and A; = 1,...,A^— 1— j. 

It remains to calculate the matrix elements Gi\fj{s,t), j = 1,. . . ,N. The inverse 
Laplace transformation of the last equation in the recursive scheme 0311) is 

GNjiw,t)= [ dt'Xit')GN-iA^-[ENit)-EMit')],t'} , j = l,...,Ar-l. (35) 

We insert herein the explicit forms of Eqs. (13 3 p and (1341) . After some algebra we finally 
obtain 



Gnn-i{w, t) = - exp {(3[w + {N - l)vt]} exp <^ -A - + (A^ - l)t, 



X 



GNj{w,t) 



X 0[-(A^ - l)vt, -{N- 2)vt; w] 



a I 

{N-2-j)\ 



N-l-j 

exp{f3[w + {N - l)vt]} J2 



1=0 



l\{N-l-j-l)\ 



X <R 



ji 



+ F,i 



w+ (N - l)vt w + {N - l)vt 
v{N~3) ' v{N-3-l) 

w + {N - l)vt 



■,w,t 



(36) 



e[-(A^ - l)vt, -{] + / - l)vt] w] + 



,t;w,t 



e[-{j + l~l)vt,-{j-l)vt;w] 



viN-j) 
Here j G {1, . . . , A^ — 2} and we have introduced the abbreviation 

Fji{a, h- w,t)= / dt' exp [-A (f, 0)] [w + {N - l)vt - (N - j - l)vt'f'^'^ . (38) 



(37) 



The main results of this Section are Eqs. (1281) and fl29l) . which give the solution of the 
Pauli rate equation ( IT2l) . and Eqs. ( I33|) -(l38l) which present the solution of Eq. ( [T71) . 
We now turn to the discussion of these results. 
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5. Discussion 

In Kittel's work |8j the equilibrium properties of the zipper are studied. The mean 
number of open links in equilibrium always increases with increasing temperature, but 
the form of this increase is different for the degeneracy factor G = 1 (cf. Fig. [3^1)) and 
for G > 1 (cf. Fig. |3|i2)). For G = 1, the mean number of open links increases smoothly 
with a concave curvature, while for G > 1, the curve exhibits a sharp increase in a narrow 
temperature interval and resembles a first-order phase transition. In the following 
discussion of representative results for the nonequilibrium dynamics and energetics, we 
use d, d/kB, and u'^ as units for energy, temperature and time, respectively. 

In the unidirectional unzipping process, the time evolution towards the completely 
unzipped state is again sensitive to the degeneracy factor G. Let us define an unfolding 
time tu by the condition that a.t t = tu the zipper has completely unfolded with a 
probability of 99,9%, i.e. 

PN{tu) = Rmitu) = I - e , (39) 

where e = 0.001. With respect to the N dependence of tu (and other quantities to be 
discussed below), we found that its behavior is similar for N = 2 and N > 2, and we 
therefore restrict the following discussion to the two-state case N = 2. Eq. ( l39l) can 
then be inverted after inserting Rjyiitu) from Eq. fl29l) in Eq. fl39l) . yielding 

tu = [l/(/3t;)] In [1 - iPv/a) In e] . (40) 

For small temperatures [large argument of the exponential and/or small prefactor 
g in Eq. (jl])] the transitions are driven predominantly by the magnitude of the energy 
gap between the closed and the opened state. For large temperatures [small argument 
of the exponential and/or large prefactor g in Eq. (jl])] by contrast, the dynamics is 
governed by the entropy difference between the closed and the opened state, and hence 
by the degeneracy factor G. The dependence of the unfolding time tu on G and T 
is plotted in FiG. [3)d) for a representative set of parameters. For any fixed nonzero 
temperature, the unfolding time decreases with increasing G, while for fixed G, the 
temperature-dependence of the unfolding time exhibits a maximum at a temperature 
T = Tniax(G), where Tmax(G') decreases with increasing G, see Fig. [3b). 

Let us now consider a certain temperature Tq and call, for this temperature, the fast- 
unzipping regime and slow-unzipping regime the ranges of G- values, where Tmax(G) < Tg 
and Tmax(G) > Tq, respectively. In these two regimes the dynamics and energetics of the 
molecular zipper exhibit a qualitatively different behavior. In particular we find (i) a 
different time-dependence of the iV-th state's occupation probability, cf. FiGS. Hb) and 
[5b), (ii) a different form of the curves describing the work needed to open the zipper, 
cf. Figs, [lb) and [5b), (iii) a different mean value of heat accepted by the zipper during 
the unzipping, cf. FiGS. [6bil) and[6bl), and iv) different values of the variances of the 
internal energy and heat during the unzipping, cf. FiGS. [6bi2) and[6b2). These features 
will be now discussed in more detail. 



13 





Figure 3. Left two panels: Mean number of open links in equilibrium for Kittel's 
molecular zipper as a function of the reservoir temperature T for A^ = 50 links, and 
(al)) G = 1 and (a2)) G = 2. Right panel b): Unfolding time iy, Eg. (gO]), as a function 
of T and G for Tq = 7.5, v = 0.25, and N = 2. Above (below) the horizontal plane 
in the graph, the opened (closed) state is energetically favored. In the base plane, the 
temperature Tmax of maximal unfolding time in dependence of the degeneracy factor 
G is shown. 



Fig. m illustrates the slow- unzipping regime. The probability PN(t) that the zipper 
has reached the opened state until time t first increases slowly. After the time d/v the 
energy of the opened state becomes lower than that of the closed one. This leads to 
more frequent transitions and accordingly pAr(t) increases more rapidly. Notice that the 
curves exhibit a change of their second derivative. The work probability density (WPD) 
during the unzipping has a maximum located inside its finite support, cf. Fig. Hb). The 
value of the work at the left (right) border of the support equals the work done on 
the zipper when it dwells during the time interval [0, tu] in the opened (closed) state. 
From the position of the WPD peak we can conclude that, for a typical trajectory of the 
stochastic process D(t), the work consists of two comparable fractions. The first (second) 
part of the work is performed while the system dwells in the closed (opened) state. At 
time tu, 0.1% of the trajectories will give molecules still residing in the zipped state. 
These trajectories contribute to the singular (5-function) components of WPDs, which 
are depicted in Fig. Hb) by the vertical arrows [331 1321 EH [35]. Note that the curves 
plotted for the temperatures T = 1 and T = 1.25, which are close to the temperature 
Traax{G) ioT G = 1, cf. FiG. |3b), become similar to the curves for these temperatures 
obtained in the fast- unzipping regime for G = 50, see Fig. |5l 

Fig. |5] illustrates the fast-unzipping regime. The probability pN(t) rapidly increases 
from the very beginning of the process, cf. FiG. |5b), i.e. the zipper opens before the 
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Figure 4. Dynamics of a) the energy levels [Eq. ([T])], b) the occupation probability 
PN{t) of the opened state [Eq. dHJ], and c) the probability density p{w,tu) of the 
work [Eq. (fT9|) ] for G == 1 and several values of the reservoir temperature T. The 
other parameters are the same as in FiG. [H The assignment of the line styles to the 
temperature given in the legend of c) applies also to parts a) and b). In a) we compare 
the dynamics of the two energy levels with the unfolding time tjj [Eq. (HH])] at different 
temperatures marked by the vertical lines. The full (dashed) line gives the energy Ei (t) 
{E2{t)) of the closed (unzipped) state. The arrows in c) represent the weights and the 
positions of the J-functions, which form the singular part of the probability density. 
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Figure 5. Same quantities as in FiG. 2] for G = 50 and otherwise the same set of 
parameters. 



15 



0.25 



-0.5^ 

"3 
> 

cij-1.25 



a1) 







— U{t) 






— W{t) 


% ^s*. 




--Q{t) 


s ^n^" 






* 

^N^ 



OJ 0.38- 



cS 0.19 
> 



a2) 



--[A[/(i)]2 

- [AW{t)f 

-[AQ{t)f 





m 0.19 4''X\ 



/« 0.13h ; 



0.06 [-• 
0, 



b2) 



1.5 



Figure 6. Upper two panels: Mean value of the internal energy [Eq. (|20p ]. work 
[Eq. dH])] and heat [Eq. ^] as a function of time for al) G = 1 and bl) G = 50. 
Lower two panels: Variances of internal energy [Eq. ((23)) ]. work [Eq. ((24)) ] and heat 
[Eq. ((25p ] as a function of time for a2) G = I and b2) G — 50. The temperature is 
T = 0.75 and the other parameters are the same as in FiG. |3| 



opened state becomes energetically preferred. The maximum of the WPD is located 
at the left border of its support, cf. Fig. [5b). This means that, for the majority of 
the trajectories, the substantial part of the work is done while the system dwells in the 
opened state. Note that the curves plotted for the temperature T = 0.25, which is close 
to the boundary temperature T^a,x{G) for G = 50, cf. Fig. [3)d), become similar to the 
curves for this temperature in the fast-unzipping regime for G = 1, see Fig. JH 

Fig. [6] illustrates the dynamics of the thermodynamic quantities (!20l) - f(25j) in the 
two unzipping regimes. For an arbitrary trajectory of D(t) which resides during the 
time interval [t',t] in the i-th state, the work performed on the system is Ei{t) — Ei{t'), 
cf. Eq. ([7]). In our model, the energies of the states decrease linearly with time and 
accordingly the mean work is a monotonically decreasing function of time, see Figs.[6|l1) 
and[6)Dl). Heat is exchanged with the reservoir when the molecule changes its state. 
It is absorbed by the molecule if the transition brings the molecular zipper to a state 
with higher energy. Since in our setting the transitions are unidirectional, the molecule 
necessarily absorbs heat up to the time tE = d/v, where the energies of the states 
become the same, cf. FiGS. H^i) andlSi). For times t > Ie, the molecule delivers heat 
to the environment. In the slow- unzipping regime we have tu > Ie- This implies that 
the mean heat first increases and then decreases, cf. Fig. [6til). By contrast, in the fast- 
unzipping regime where tu < tE, the mean heat monotonically increases, cf. Fig. [Gbl). 
Finally, due to the transitions to the state with higher energy at the very beginning of 
the process, cf. FlGS.|Hand[5l the mean internal energy ( |20l) develops a single maximum. 
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The variances of U(t), W(t) and Q(t), cf. Eqs. (^-(^, are plotted in FiGS. [Si2) 
and Figs. [6)32). If the variance of the internal energy approaches zero, the variance 
of the work becomes equal to that of the heat. All variances approach constant values 
at large times. In fact, for t > tu, almost all trajectories have brought the molecule 
into the opened state. As a consequence, the increments to the work, heat and internal 
energy are nearly constant during the time interval [t(7,t]. Therefore the form of their 
probability densities does not change, the curves just move along the energy axis. 

The probability density for the internal energy has no continuous component. In 
general it consists of A^ ^-functions located at the energies of the individual states. 
The corresponding weights are given by the occupation probabilities of the states. In 
the slow-unzipping regime, [Af/(t)]^ vanishes at time t = and at time t^ = d/v, 
when the state energies are equal, and it becomes very close to zero for t > tu- In 
between these time instants, the variance develops a maximum, cf. Fig. [6ti2). In the 
fast-unzipping regime, the function [Af/(t)]^ displays just one maximum, cf. Fig. [6]d2), 
because tu <tE- 

In both the fast and slow unzipping regime, the variance of the work monotonically 
increases and approaches a constant for large times, cf. Fig. [6|l2) and Fig. |6b2). The 
absolute value of the first derivate of the variance [/S.Qit)]^ is given by the product 
of the two quantities, the energy difference between the closed and the opened states 
at time t and the probability that the zipper opens at the time t, which were shown 
in Figs. [Hand [51 In the slow-unzipping regime, the majority of the transitions occurs 
during the time interval [tg, tu\- During this time interval, the energy difference between 
the states monotonically increases. In the fast-unzipping regime, by contrast, nearly all 
transitions take place before time t£;/2. The sooner the transition, the larger is the 
amount of transferred heat. As a result, at small times, the heat probability density has 
one peak at a large value, coming from the trajectories with a transition to the opened 
state, and another peak at zero heat exchange, originating from trajectories without a 
transition. With increasing time the amount of the trajectories with a transition to the 
opened state increases rapidly, and the peak close to zero heat moves towards the peak 
at a large heat value. This explains the behavior of the heat variances in FiGS. [6^2) 
and[6)o2). 

Finally, one may ask how our results are affected if the variation ( "static disorder" ) 
in base pairing energies is included in the modeling. To this end we have performed 
Monte Carlo simulations |18] of the stochastic process for a molecular zipper with 
A^ = 10 states, where the initial energies i?fc(0), /c = 1, . . . , A^, in Eq. ([T]) are given by 
Ek{0) = (A; — l)(A + ?7fc), corresponding to different losses of energies due to variations in 
base pair bondings. The r)k were chosen as random numbers from a box distribution in 
the interval [—A/3, A/3]. Results from these simulations for a number of realizations of 
this disorder were compared to the predictions of the analytical theory for the "ordered 
case", where -E'fc(O) = {k — 1)A. We found that, for typical realizations of sets of rjk, 
the shapes of the work distributions are very similar, while the peak positions and peak 
heights are shifted slightly. Also the probability distributions PN{t) for the zipper to 
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fully open until time t are, for these sets of 77^, very similar in shape. Analogous to 
the peak positions of the work distributions, the onset of opening shifts slightly from 
realization to realization. The shifts of the peaks and of the onset of the opening are 
controlled by the largest base pair bonding (largest r^^), which governs the unfolding 
time. 
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